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Nuclear reactors provide intense sources of electron antineutrinos, characterized by few-MeV 
energy E and unoscillated spectral shape &(E). High-statistics observations of reactor neutrino os¬ 
cillations over medium-baseline distances L ~ 0(50) km would provide unprecedented opportunities 
to probe both the long-wavelength mass-mixing parameters ( 5m 2 and 0 12 ) and the short-wavelength 
ones (A m 2 e and # 13 ), together with the subtle interference effects associated with the neutrino mass 
hierarchy (either normal or inverted). In a given experimental setting—here taken as in the JUNO 
project for definiteness—the achievable hierarchy sensitivity and parameter accuracy depend not 
only on the accumulated statistics but also on systematic uncertainties, which include (but are not 
limited to) the mass-mixing priors and the normalizations of signals and backgrounds. We examine, 
in addition, the effect of introducing smooth deformations of the detector energy scale, E —> E'(E), 
and of the reactor flux shape, &{E) —> &(E), within reasonable error bands inspired by state-of- 
the-art estimates. It turns out that energy-scale and flux-shape systematics can noticeably affect 
the performance of a JUNO-like experiment, both on the hierarchy discrimination and on preci¬ 
sion oscillation physics. It is shown that a significant reduction of the assumed energy-scale and 
flux-shape uncertainties (by, say, a factor of 2) would be highly beneficial to the physics program of 
medium-baseline reactor projects. Our results also shed some light on the role of the inverse-beta 
decay threshold, of geoneutrino backgrounds, and of matter effects in the analysis of future reactor 
oscillation data. 


I. INTRODUCTION 

The currently established neutrino oscillation phenomenology can be interpreted in a three-generation framework, 
whose relevant parameters are three mixing angles ( 812 , $ 13 , $ 23 ), one possible CP-violating phase S, and two neutrino 
squared mass differences { 8 m 2 , ±Am 2 ), where the latter sign distinguishes the cases of normal hierarchy (NH, +) and 
inverted hierarchy (IH, —) for the neutrino mass spectrum jj]. In reactor neutrino experiments, for a given neutrino 
energy E and baseline L , the V e disappearance probability depends only on a subset of parameters, 

P{v e ->• V e ) = P ee ( 0 i 3 , $ 12 , 8 m 2 , ±Ato 2 ) , (1) 

and is generally not invariant under a swap of the Am 2 sign, thus providing some sensitivity to the mass hierarchy 
as noted in (2j. Hierarchy effects vanish in the limits of short baselines [L < 0{ 1) km] and of long baselines [L > 
0(1O 2 ) km], where the P ee arguments reduce to ($ 13 , |Am 2 |) and ($ 13 , $ 12 , 8 m 2 ), respectively. However, at medium 
baselines L (few tens of km), the hierarchy sensitivity can be recovered in high-statistics reactor experiments, provided 
that interference effects between short- and long-wavelength oscillations are resolved, as originally proposed in [3]. 

Such a possibility is being seriously investigated experimentally, as one of the main physics goals of the Jiangmen 
Underground Neutrino Observatory in China [ 3 ] (JUNO, in construction), and of the project RENO-50 proposed in 
Korea [5j]. A number of papers have examined the stringent conditions needed to discriminate the hierarchy in this 
class of reactor experiments, including high energy resolution, small dispersion of multiple core baselines, good control 
of energy-scale nonlinearities, and reduction of systematics related to oscillation parameters and normalizations. 
Satisfying these conditions would also provide more precise measurements of the parameters ($ 12 , 8 m 2 ) 10 and of 
Am 2 . We refer the reader to Si for recent reviews of this field of research and of the related bibliography. 

In this context, it has been realized that the shape of the (unoscillated) reactor neutrino flux &{E) may not be 
known with the desirable accuracy, as recently demonstrated by the unexpected spectral features consistently found 
in the current short-baseline experiments RENO 0 , Double Chooz [ll| and Daya Bay 0 . In particular, an event 
excess (sometimes dubbed as a “bump” or “shoulder” in the spectrum) clearly emerges around E ~ 5-7 MeV, with 
respect to widely adopted theoretical predictions [3 EH- Very recent calculations provide possible (although still 
partial) interpretations of the subtle nuclear effects which may be responsible for these (and possibly other) unexpected 
spectral features p~5Ml~7l ] • but a thorough understanding of their origin has not yet been achieved. 

In current short-baseline experiments, endowed with near and far detectors, poorly understood spectral features 
largely cancel in near-to-far flux ratios, and do not significantly affect the measurement of the dominant parameters 
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(013, |Am 2 |) [IS] - However, in the absence of a near detector, a medium-baseline experiment such as JUNO must 
rely on absolute estimates of the flux spectrum $(U) and of its associated uncertainties, which may affect both the 
sensitivity to the hierarchy and the accuracy of the ($ 12 , 5m 2 , |Am 2 |) measurements. 

A related and subtle issue concerns energy-scale variations E —► E'(E) [Ij| which, in principle, may be nonlinearly 
engineered to produce a hierarchy ambiguity [20j . Although the energy-scale issue can be kept under control by 
calibration constraints at subpercent level [2jJ, specific combinations of energy variations E —> E'(E) and spectral 
deviations $(£!) — > $'(£') might represent a more subtle threat to the hierarchy discrimination f 22j . In general, deli¬ 
cate statistical aspects—related to the treatment of admissible spectral deformations—are now emerging in neutrino 
oscillation searches, mirroring the evolution of other fields of physics from the discovery phase to the precision era, as 
remarked in [23[ . 

Within this scenario, and building upon our previous work [ 2 ^ ]. we present herein a systematic analysis of the 
combined effects of energy-scale nonlinearities E —> E'(E) and flux-shape uncertainties $(£?) —> &'(E), assuming a 
reference JUNO-likc medium-baseline reactor neutrino setting. The structure of our paper is as follows: In Sec. II 
we describe the adopted notation and methodology. In Sec. Ill we discuss the effects of energy-scale and flux-shape 
uncertainties on the hierarchy discrimination, while in Sec. IV we show their impact on precision measurements of 
the oscillation parameters. In Sec. V we repeat the analysis in a prospective scenario where the energy-scale and flux- 
shape uncertainties are reduced by a factor of two. We also comment on the role of the inverse-beta decay threshold, 
geoneutrino backgrounds, and matter effects. A brief summary of the results, and perspectives for further work, is 
presented in Sec. VI. 


II. NOTATION AND METHODOLOGY 


We generally adopt the same notation and inputs as in 1221 (to which the reader is referred for details), and comment 
below only about those aspects that are new, modified, or relevant for the present analysis. In particular, we discuss 
the statistical approach used to characterize and include energy-scale and flux-shape uncertainties. We remark that 
the adopted experimental setup [ 
the JUNO Collaboration 01. 


as taken from [21], is basically the same as reported in a recent publication from 


A. Neutrino oscillation parameters and priors 


Neutrino oscillation probabilities can be expressed in terms of neutrino squared mass differences (Am 2 , = m 2 — ?n 2 ) 
and trigonometric functions of the mixing angles 0^ (e.g., Sj 2 = sin 0, j , c,; ? = cos 0, j ). In medium-baseline reactor 
experiments, there are four relevant parameters: s 2 2 , sf 3 , 5m 2 = Am 2 !, and A m 2 e defined as [mI - IaI 


Am 2 e = Am 2 ± ~(c? 2 - s 2 12 )5m 2 , (2) 

where the upper (lower) sign refers to NH (IH), while Am 2 is defined as nnm 

Am 2 = — | Am. 2 ! + Am 32 1 > 0 . (3) 


An accurate analytical expression for the oscillation probability P ee , including effects due to propagation in constant- 
density matter and to multiple reactor cores, can be found in [ 22 } [see Eqs. (58) and (59) therein]. 

We assume the following reasonable priors (central values and ±ler errors) for the above parameters, at the start 
of a JUNO-likc experiment: 


s \ 2 = (3.08 ±0.17) x 10" 1 , (4) 

Sm 2 = (7.54 ±0.20) x 10" 5 eV 2 , (5) 

s 2 13 = (2.20 ±0.08) x 10" 2 , (6) 

Am 2 e = (2.40 ± 0.05) x 10" 3 eV 2 . (7) 


The (s 2 2 , 5m 2 ) priors—unlikely to change significantly in the near future—are taken from the global fit in (27} . with 
errors defined as 1/6 of the ±3cr range. The sf 3 central value is a bit lower than in [27], as suggested by recent reactor 
results ! tutu, and is also endowed with a smaller ±ler error, representative of the final accuracy expected in Daya 
Finally, the A m 2 e central value is also in ballpark of the current global fits [l^|]^,[3Cj, but with a somewhat 


Bay 

smaller fractional error than in [27| (2.0% instead of 2.6%) as it can be expected from near-future improvements in 
short-baseline reactor [28} and long-baseline accelerator experiments [3l} . 
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B. Reactor and geoneutrino spectra, energy resolution and thresholds 


Concerning the reactor neutrino fluxes (from both medium-baseline and far sources), we use the same average 
fuel component and overall normalization as in [22} . but we alter the energy profile to include the newly discovered 
spectral feature at E ~ 5-7 MeV jTol - lTl ]. In particular, we multiply the unoscillated reactor spectrum in [22} by a 
smoothed version of the bin-to-bin ratio (Daya Bay)/(Huber ± Mueller) reported in j28}, which effectively accounts 
for the spectral bump feature as observed in Daya Bay [H}. 

The normalization of the U and Th geoneutrino background fluxes in JUNO is slightly increased with respect to 
(22} by 7% and 9%, respectively, in order to match the most recent estimates in [!2[. We also update the associated 
errors, by conservatively taking the largest of the asymmetric lcr uncertainties from 32], namely, ±20% for the U flux 
and ±27% for the Th flux. 

Finally, with respect to j22s], we slightly update the energy resolution width a e as reported in [28} . 


cr e (E e ) 
E e ± m e 


2.57 x 10" 2 
\f(E e ± m e )/MeV 


±0.18 x 10“ 2 , 


( 8 ) 


where E e ± m e is the total (true) visible energy of the inverse beta decay (IBD) event, as a sum of the total positron 
energy E e and electron mass m e . We remind the reader of the following 1221] : (1) the finite resolution width smears 
the observed visible energy E v - ls around its true value E e ± m e ; (2) the neutrino energy E has an IBD kinematical 
threshold E > 1.806 MeV; (3) the parent neutrino energy E and the observed visible energy of the event E v ; s are 
approximately related by 


E ~ E vis ± 0.78 MeV , (9) 

up to nucleon recoil and energy resolution effects (which we accurately include in the calculations of energy spectra 
[23]); and (4) the visible energy threshold is thus E v - ls > 1 MeV. 


C. Error bands for energy-scale and flux-shape deformations 


In first approximation, we assume a JUNO energy-scale uncertainty comparable to the Daya Bay one. For current 
Daya Bay data, the la error band of admissible deviations in the reconstructed/true visible e nergy ratio has been 
shown in |l8[ and (with slightly smaller width) in (33[, H|J. We have translated the bands in [33, Hi into relative 
deviations E'/E for the neutrino energy via Eq. ([9]). Asymmetric lcr uncertainties have been symmetrized to the 
largest between ±lcr and —lcr. Figure 1 (top panel) shows, in color, the resulting energy-scale error band (at ±lcr 
in E'/E), as a function of the parent neutrino energy E. Besides this “default” band, we shall also consider a 
more optimistic case with “halved” errors (dot-dashed lines in the top panel of Fig. 1), in view of dedicated energy 
calibration campaigns expected in JUNO. 




FIG. 1: Default zblcr error bands assumed for energy-scale deviations E'/E (top panel) and flux-shape variations (bottom 

panel), in terms of the neutrino energy E. Bands with halved errors are also shown (dot-dashed lines in both panels). 
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Concerning the flux-shape uncertainties of the unoscillated reactor spectrum Q(E), we assume that (5>'(E)/$(E) 
deviations are constrained by the ± 1 <t error bands estimated in [ 15j. We have smoothed out and symmetrized the 
bands in a as reported in Fig. 1 (bottom panel) in terms of the neutrino energy E. Since the issue of reactor 
spectral shapes is still highly debated, the error band should be taken as merely indicative of the current 

level of theoretical uncertainties. The high statistics accumulated in the present generation of short-baseline reactor 
experiments will certainly help to constrain any model of reactor spectra and, indeed, the current size of systematic 
shape uncertainties estimated in Daya Bay [l 8 | already seems to be a factor of 2 smaller than in [l5j], although a 
detailed assessment has not yet been published. For this reason, also in the analysis for flux-shape uncertainties, we 
shall consider the more optimistic case of halved theoretical errors (dot-dashed lines in the bottom panel of Fig. 1). 

A final remark is in order. In the absence of a detailed characterization of the error bands in Fig. 1, we simply assume 
that they scale linearly with ncr. We also neglect, for lack of published information, possible error correlations at 
different energies. Although some correlations are known to exist, as a result of underlying models for both the energy 
scale nonlinearities pl| HU, HI] and the reactor spectra fl3l - [l6| , their impact should not be overemphasized at this 
stage. Indeed, the recently observed, localized “bump” feature lar gely exceeds the estimated errors and covariances 
that where thought to characterize the spectra a few years ago [13|, Qjj. In this sense, neglecting possible covariances 
in Fig. 1 should lead to conservative results. A more refined analysis will be possible when such error bands will be 
determined more precisely, and endowed with point-to-point correlation functions. 


D. Statistical approach 


As in [22|, we calculate the “true” event spectrum S*(E V i s ) by assuming the central values of the oscillation 
parameters reported in Sec. II A, for either NH or IH. Such a spectrum S* represents the “experimental data”, to 
be compared with a family of spectra S(E v - ls ), obtained by varying the continuous parameters ( dm , Ara^, $ 12 , $ 13 ), 
in either the same or the opposite hierarchy. The comparison is performed in terms of a y 2 function that contains 
statistical, parametric, and systematic components (see also H), 


2 2 , 2,2 

X Xstat A Xpar ~b X; 


sys 


( 10 ) 


The Xstat term (which embeds statistical fluctuations) is the same as in 22]. The y 2 term (which embeds penalties 


for the oscillation parameters) is also unchanged, apart from the numerical priors, here taken from Sec. II A above. 
The Xsys term contains, as in |22], two penalties related to the geoneutrino flux normalizations, 


Xgeo 


i=u,Th 


fj -1 


e Xs 


(ii) 


where we take sy = 0.20 and syh = 0.27, as described in Sec. II B. However, while in 22] the y 2 ys term was completed 
by just another penalty for the overall reactor flux normalization, now it must be supplemented by appropriate 
penalties for energy-scale and flux-shape deformations. 

To this purpose, we consider smooth deformations of the energy scale E E'(E) (that we assume to act upon 
the “experimental spectrum” S*) and of the flux shape <1>(.E) —> &(E) (that we assume to act upon the “theoretical 
spectrum” S), in terms of generic polynomials in E (in MeV), 


k 

i + J2^ Ei = 1 + M#) , (12) 

2=0 
h 

l + J2PjE j = 1 + 6*(E) , (13) 

i=o 

with h and k increasing until stable results are reached [36|. Note that the trivial cases h = 0 and k = 0 correspond, 
respectively, to an overall renormalization of the energy scale [E' = (1 + ao)E] and of the reactor spectrum [$' = 

(1 + Ao)$]- 

With reference to Fig. 1, let us denote the boundaries of the lcr error bands in Fig. 1 as 1 ± Se{E ) for the upper 
panel, and as 1 ± S$(E) for the lower panel. Then we define two new systematic penalties, in terms of the largest 
relative deviation associated with each polynomial: 


E' 

~E 

&(E) 

${E) 


Xe = 


8e(E) 

Se(E) 


max 

E 


(14) 
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x! = max 

h/ 

In other words, if the polynomial function 5e{E) “touches” the na error band boundary n x Se{E ), its contribution 
to the Xsys term is assumed to be n 2 , and similarly for 5&(E) and S^(E). Equivalently, the ±ler bands in Fig. 1 
are assumed to be the envelope of all possible systematic deviations at the ler level, and similarly for na. Such a 
y 2 characterization of energy-scale and spectral-shape errors is both intuitive and conservative, as appropriate to an 
exploratory analysis. As previously remarked, more refined definitions of y 2 ys will be possible in the future, in terms 
of energy-dependent cross correlations. 

We have found that our results, to be discussed in the next sections, become numerically stable already for fifth-order 
polynomials, which are taken as a default choice for all the following figures. Therefore, in general, the y 2 minimization 
requires scanning a 18-dimensional parameter space, including four oscillation parameters (s 2 2 , sf 3 , 5m 2 , Am 2 e ), two 
geoneutrino flux normalizations (/u, /rh), and twelve polynomial coefficients (op,..., as) and (/?o,..., /3s). [We have 
also cross-checked the numerical results by using different and independent minimization methods.] 

For the sake of the discussion, we shall also consider cases with reduced dimensionality, as obtained by setting to 
zero the coefficients cn or /3j. However, in our analysis, the specific coefficient (3q is never zeroed a priori, since it 
parametrizes a floating normalization for the reactor flux, $ —> <f>(l + /3o)- In particular, we shall consider the following 
cases, in order of increasing number of free parameters: 

• oscillation + normalizations: {s\ 2 , s 2 3 , 5m 2 , A m 2 e ) + (/u, /Th) + (/?o) = 7 parameters; 

• osc. + norm. + energy scale: as above + (ao,..., as) = 13 parameters; 

• osc. + norm. + energy scale + flux shape: as above + (/3i,..., /3s) = 18 parameters. 

In the first two cases, from the definition of y|,, the lcr error associated with /3q coincides to the smallest error band 
width in Fig. 1 (bottom panel), i.e. to ~ 2.3%, which is a typical value for the reactor flux normalization uncertainty. 


6*(E) 


S*(E) 


(15) 


III. DEFAULT ENERGY-SCALE AND FLUX-SHAPE ERRORS: HIERARCHY TESTS 


In this section we study the hierarchy sensitivity of the reference JUNO experiment, under the effects of increasingly 
large sets of systematic errors. Default lcr errors are assumed for energy-scale and flux-shape deviations (see Fig. 1). 
The various effects are first shown graphically and then quantified in terms of the variable 37] N a = \J Ay 2 


Figure 2 (upper panel) shows the absolute NH and IH event spectra S NH and S'™, respectively, as obtained by simply 
swapping the hierarchy, without any change in the central value of the oscillation parameters or other systematics 
(case of “no uncertainties”). For definiteness, the spectra correspond to 5 years of data taking. The lower panel 
of Fig. 2 shows the IH/NH spectral ratio, with characteristic wiggles due to the mismatch between the oscillation 
peak positions in the two hierarchies. The amplitude of the wiggles reaches ~ 4% for E v - ls ~ 2-3 MeV. In this ideal 
situation, the NH and IH spectra would be distinguishable “by eye.” 

In Fig. 3, the NH spectrum is taken as input data (S* = S' 1 " 111 ), while the IH spectrum is fitted (S = S IH ), with 
allowance for oscillation parameter uncertainties and normalization systematics. With respect to Fig. 2, the mismatch 
between NH and IH spectra appears to be reduced to < 3% in the low-energy range, although it slightly increases 
(up to ~ 1.5%) at high energy. 

Figure 4 is similar to Fig. 3, but with energy-scale uncertainties included in the fit. In this case, the NH and IH 
spectra are barely distinguishable by eye, their relative mismatch being lower than ~ 2% at any energy. This trend is 
even more pronounced in Fig. 5, where flux-shape uncertainties have been included in the fit: the NH and IH spectra 
appear to be almost indistinguishable, except for percent-level differences in the oscillation peaks around 2 MeV. 

In Figs. 2-5, a sharp rise of the event spectra is evident at the IBD threshold. The steep derivative guarantees that 
energy-scale and flux-shape deviations cannot be large at threshold, otherwise the spectral mismatch would locally 
“explode,” with a significant y 2 increase. In this sense, the IBD threshold acts as a “self-calibrating point”: the true 
(NH) spectrum and the test (IH) spectrum must almost coincide at such kinematical threshold, and systematic errors 
in the fit cannot alter this requirement. [At most, the small residual mismatch around threshold may lead to locally 
fuzzy wiggles in the spectral ratio S*/S, as partly shown in Figs. 4 and 5.] 

The above expectations are confirmed by an analysis of the best-fit energy profiles (fifth-degree polynomials) for 
the energy-scale and flux-shape deviations, corresponding to the cases shown in Figs. 3, 4 and 5. Figure 6 shows such 
profiles (solid curves) superimposed to the default error bands (in color) for E'/E (top panels) and $'/$ (bottom 
panels). The leftmost panels of Fig. 6 correspond to the fit in Fig. 3, which includes only oscillation and normalization 
errors. In this case, E'/E = 1 by construction (no energy-scale error), while $'/4> = 1 + /3q can float (to account for 
the flux normalization), but happens to have a best-fit value very close to unity. 
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FIG. 2: Comparison of event spectra in NH and IH, as obtained 
for fixed oscillation parameters and no systematic errors. Top: 

A Konluto ono/’tro f/~»v '7 1 — ^ irr Rn+tom ■ Qnocf rol rotin 



FIG. 3: Comparison of event spectra in NH (true, S*) and IH 
(fitted, S), including oscillation and normalization uncertainties. 




FIG. 4: As in Fig. 3, but including energy-scale systematics. 


FIG. 5: As in Fig. 4, but including flux-shape systematics. 


The middle panels in Fig. 6 correspond to the fit in Fig. 4, which also includes energy-scale systematics. In this case, 
one can observe a slight offset of the overall ratio <h'/4> = 1 + /?o, and a peculiar pattern for the best-fit E'/E profile. 
The function E'/E is close to unity at the IBD threshold (E ~ 1.8 MeV, equivalent to E v ; s ~ 1 MeV), according 
to expectations. Then the function rises up by ~ +0.6%, changes sign at E ~ 4 MeV, decreases by ~ —0.6%, and 
approaches unity at high energy. The sign-flip behavior is vaguely reminiscent of the “engineered” E'/E profile that 
would realign, by construction, the IH and NH oscillation phases and peaks [20] while remaining unity at the IBD 
threshold (see Fig. 17 in [23). However, the “fitted” E'/E profile in Fig. 6 does not need to reach large deviations 
at high E as the engineered one f+j . since the high-energy tail of the spectrum contributes marginally to the y 2 
function. Finally, the rightmost panels in Fig. 6 correspond to the complete fit in Fig. 5, which includes also flux- 
shape systematics. The function E'/E is qualitatively similar to the middle panel, but with reduced deviations in 
the high-energy part of the spectrum. The best-fit function $'/$ shows sign-changing deviations at the few-percent 
level, well within the ±lcr (colored) error band. In conclusion, admissible systematic deformations of the energy scale 
and of the flux shape, added to the usual oscillation parameter and normalization uncertainties, may bring the “true” 
and “wrong” event spectra as close to each other as shown in Fig. 5. 
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osc. + norm. + energy scale + flux shape 



FIG. 6: Energy profile of best-fit deviations E'/E (top panels) and 'F/F (bottom panels), for different sets of systematic uncertainties. 


Figure 7 shows the statistical significance of the wrong hierarchy rejection, for the case of true normal hierarchy, in 
terms of N a = yj Ay; 2 as a function of live time T. The abscissa scales as y/T, thus showing at a glance any deviation 
from the ideal “linear” case of purely statistical errors (N a oc y/T). In the fit including only oscillation parameter and 
normalization uncertainties, N a grows steadily and almost linearly in y/T along ten years of data taking. However, 
the inclusion of energy-scale uncertainties provides some bending of the linear rise, with a noticeable but not dramatic 
decrease of the statistical significance. In particular, it appears that a 3er rejection is achievable after about six years 
of data taking. We agree with [21] that energy-scale uncertainties, by themselves, do not represent a showstopper for 
JUNO-like experiments. However, Fig. 7 shows that the combination of energy-scale and flux-shape systematics can 
be quite sizable: the solid curve for N a grows much more slowly than y/T, and remains below 3 a even after ten years 
of data taking. Figure 8 shows a very similar behavior, but assuming the IH as true. Figures 7 and 8 represent one 


Na 



Na 



FIG. 7: Case of true NH: Statistical significance of the IH rejec¬ 
tion as a function of the detector live time T, as derived from fits 
including different sets of systematics. Note that the abscissa 
scales as y/T. The horizontal 3cr line is shown to guide the eye. 


FIG. 8: As in Fig. 7, but for true IH and rejection of NH, 

























FIG. 9: Mass-mixing parameters (5m 2 , s 2 2 ): l<r contours for 
true NH and T = 5 y, as derived from fits including different 
systematic uncertainties. The arrows indicate the best-fit dis¬ 
placement in the cases of wrong hierarchy (green) and no matter 
effects (magenta). 



FIG. 10: As in Fig. 9, but for the (Am 2 e , s 2 3 ) parameters. 


IV. DEFAULT ENERGY-SCALE AND FLUX-SHAPE ERRORS: PRECISION PHYSICS 

In this section, unless otherwise noticed, the true hierarchy is assumed to be known, and the discussion is focused on 
the precision physics program, namely, on the accuracy expected at Tier on the oscillation parameters and geoneutrino 
normalizations. For definiteness, the accumulated statistics refers to T = 5 yr. 

Figure 9 shows the ler contours in the plane charted by the mass-mixing parameters (Sm 2 , sf 2 ), assuming true NH. 
In the case with only oscillation and normalization errors, the accuracy of both Sm 2 and s\ 2 is more than an order 
of magnitude better than the prior errors assumed in Sec. II A. The accuracy is slightly worse in the presence of 
energy-scale systematics, but is significantly degraded (by almost a factor of three) in combination with flux-shape 
systematics. This is not surprising, since the (Sm 2 , s 2 2 ) parameters govern the long-wavelength oscillation pattern, 
which is sensitive to smooth deformations of the spectral shape. Note that, in all the above cases of Fig. 9, the best-fit 
coordinates coincide—by construction- with the central values of (Sm 2 , s 32 ) assumed as priors. However, the best-fit 
point would be significantly displaced if the inverse hierarchy were mistakenly assumed as “true,” as indicated by the 
green arrow in Fig. 9 (color online). A different but still sizable displacement (indicated by the magenta arrow) would 
also be induced by discarding matter effects in the fit. For simplicity, the matter density has been assumed herein to 
be constant [ 22 j but, in the future, one should characterize more precisely the density profile from the geophysical and 
geochemical viewpoint. Summarizing, Fig. 9 shows that the ler accuracy of the (Sm 2 , Sl2 ) measurements in a JUNO- 
like experiment can be significantly degraded by the combined effect of energy-scale and flux-shape uncertainties, and 
that their central values may be biased by ambiguities in the hierarchy, as well as by “vacuum” approximations in 
the oscillation probability. 

Figure 10 shows the ler contours in the plane charted by the mass-mixing parameters (Am 2 e , s 2 3 ), assuming true 
NH. The accuracy on s 33 is essentially constant and almost equal to the prior assignment in Sec. II A, implying that 
a JUNO-like experiment cannot really improve the input 6 13 data from current short-baseline reactor experiment. In 
the case with only oscillation and normalization errors, the accuracy of A m 2 e is more than an order of magnitude 
better than the prior error assumed in Sec. II A. However, the accuracy is degraded (by a factor of two) by energy- 
scale uncertainties and (by a factor of three) by adding flux-shape systematics. In fact these systematics, as shown 
in Sec. Ill, may slightly alter the pattern of short-wavelength oscillations and may thus affect the measurement of its 
peak frequency, governed by Am 2 e . Moreover, the central value of A m 2 e (but not of # 13 ) would be strongly biased 
in the case of “wrong” hierarchy (green arrow), roughly by (c 2 2 — s 2 2 )Sm 2 as expected from Eq. J2]). The vacuum 
approximation bias (magenta arrow) is instead insignificant, since the parameters (Am 2 e , s 33 ) are basically unaffected 
by matter effects [22|. Results similar to Figs. 9 and 10 also hold in the case of true IH (not shown). 
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FIG. 11: As in Fig. 9, but for (Th, U) geo-// normalizations. 


Finally, Fig. 11 shows the ler contours in the plane charted by the Th and U geoneutrino flux normalizations 
(/Th, /u)- The fit constrains these normalizations within ler errors which are smaller (by about 30%) than their prior 
values as defined in Sec. II D, and are quite insensitive to different sources of systematic errors and biases. Therefore, 
prospective geoneutrino results in JUNO will help to constrain better the current geophysical and geochemical models 
for the radiogenic element abundances, independently of systematic details. 


V. RESULTS FOR HALVED ENERGY-SCALE AND FLUX-SHAPE ERRORS 

The reference error bands in Fig. 1 are representative of present state-of-the-art systematic uncertainties on energy- 
scale and flux-shape deformations of typical reactor event spectra. When one or more JUNO-like experiments will 
be operative, it is conceivable that the detector energy scale will be subject to dedicated calibration campaigns, and 
that the reactors spectral profiles will be better understood, on the basis of the high-statistics data sets collected by 
current-generation short-baseline experiments. Therefore, it makes sense to repeat the analyses in Secs. Ill and IV 
in the hypothesis of smaller (for definiteness, halved) error bands in Fig. 1, while all other priors are assumed to be 
unchanged. The results are discussed below. 



FIG. 12: As in Fig. 5, but with halved energy-scale and flux-shape errors. 
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FIG. 13: As in Fig. 6, but with halved energy-scale and flux-shape uncertainties. 

A. Hierarchy tests 

Figure 12 is analogous to Fig. 5, but with halved energy-scale and flux-shape uncertainties. It can be seen that the 
NH and IH spectra cannot be brought as close to each other as in Fig. 5, and that the residual spectral differences 
are noticeably larger, as a result of the systematic error reduction. 

Figure 13 shows the effect of halving errors on the energy profiles of the energy-scale and flux-shape deformations, 
corresponding to the best-fit IH spectrum in Fig. 12. From the comparison of Fig. 6 and Fig. 13, it appears that the 
qualitative behavior of the profiles is similar for either default or halved errors but, in the latter case, the amplitude 
is suppressed by roughly a factor of two. 

From Figs. 12 and 13, one expects a non-negligible improvement in testing the wrong hierarchy versus the true 
one. Indeed, Figs. 14 and 15 show the statistical significance of the wrong-hierarchy rejection, in the cases with true 
NH and true IH, respectively: in both cases, a 3er rejection level appears to be reachable in about 6 years of data 
taking, consistently with the expected goal of a JUNO-like experiment f2lf. By comparing Figs. 14 and 15 with the 
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FIG. 14: As in Fig. 7 (true NH), but with halved energy-scale FIG. 15: As in Fig. 14, but for true IH. 

and flux-shape uncertainties. 
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FIG. 16: As in Fig. 9, but with halved energy-scale and flux- FIG. 17: As in Fig. 10, but with halved energy-scale and flux- 
shape uncertainties. shape uncertainties. 


current estimates, may prevent (in combination) an effective hierarchy discrimination; (3) a future error reduction by 
a factor of two may lead to a > 3<r rejection of the wrong hierarchy, with a reasonable detector exposure (T ~ 6 y); 
(4) a significance > 4 <j seems to be out of reach for a one-decade exposure, unless all systematics are further reduced. 

B. Precision physics 

Halving the energy-scale and flux-shape uncertainties has also a significant impact on the precision program in a 
JUNO-like experiment. Figures 16 and 17 (to be compared with the analogous Figs. 9 and 10, respectively), show that 
such uncertainties, with respect to the case with only oscillation and normalization errors, degrade the final accuracy 
on the ( 5m , sf 2 , A rn ^ e ) parameters by a factor of two or less, while the sj 3 accuracy remains always close to its prior 
assignment. The constraints on the Th and U geoneutrino flux normalizations (not shown) are basically the same as 
in Fig. 11, since these spectral components are quite insensitive to details of the energy scale and flux shape. 

Table I summarizes the fit results for the oscillation and the geoneutrino parameters, in terms of (symmetrized) 
lcr errors, to be compared with their prior ±1 ct ranges. The results refer to NH and to increasing sets of systematic 
uncertainties, including energy-scale and flux-shape errors at default values (see Figs. 9 and 10) and halved values (see 
Figs. 16 and 17). In conclusion, reducing the error bands in Fig. 1 represents a major requirement to fully exploit the 
physics potential of a JUNO-like experiment, both for discriminating the neutrino mass hierarchy and for measuring 
more precisely the ( 8m , s 32 , Am^) oscillation parameters and the geoneutrino fluxes. 


TABLE I: Precision physics in a JUNO-like experiment, assuming known normal hierarchy. 1st and 2nd column: oscillation or geoneutrino 
parameter, together with the assumed prior value and ±1 a error. 3rd column: 1 a error from the fit to prospective 5-year data, including only 
oscillation and normalization uncertainties. 4th and 5th column: la error from the fit, including also energy-scale and flux-shape uncertainties 
with default error bands. 6th and 7th column: as in the previous two columns, but with halved error bands. Similar results are obtained for the 
case of known inverted hierarchy (not shown). See the text for details. 


Parameter 

Prior ±1(7 

Osc. + norm. 

fit error 

+ Energy scale 
(default) 

+ Flux shape 
(default) 

+ Energy scale 
(halved) 

+ Flux shape 
(halved) 

s^/10" 1 

3.08 ± 0.17 

0.015 

0.021 

0.040 

0.017 

0.026 

(5m 2 /10 -5 eV 2 

7.54 ±0.20 

0.016 

0.017 

0.038 

0.016 

0.029 

s i 3 /io ~ 2 

2.20 ± 0.08 

0.073 

0.073 

0.074 

0.074 

0.074 

Am 2 e /10- 3 eV 2 

2.40 ± 0.05 

0.0036 

0.0074 

0.011 

0.0064 

0.0074 

/Th 

1.00 ± 0.27 

0.20 

0.21 

0.21 

0.21 

0.21 

/u 

1.00 ± 0.20 

0.14 

0.14 

0.14 

0.14 

0.14 
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VI. SUMMARY AND CONCLUSIONS 

Medium baseline, high-statistics reactor neutrino projects such as JUNO (in construction) and RENO-50 (proposed) 
can pursue an important research program in neutrino physics, including the determination of the unknown mass 
hierarchy, precision measurements of some known oscillation parameters, and improved constraints on geoneutrino 
fluxes. In this context, building upon our previous work [22j , we have examined in detail the effects of nonlinear 
variations of the energy scale, E — > E'(E), and of the unoscillated reactor neutrino flux shape, $(£?) — > $>'(E), in 
addition to the usual prior uncertainties associated with oscillation and normalization parameters. For definiteness, we 
have performed our analysis in a JUNO-like configuration, assuming energy-scale and flux-shape error bands anchored 
to state-of-the-art estimates (default case), as well as for error bands reduced by a factor of two (halved case), as 
shown in Fig. 1. 

It turns out that such systematics can noticeably affect the performance of the experiment, and that their reduction 
is mandatory in order to achieve statistically significant results, both for hierarchy discrimination and for precision 
physics. In particular, a > 3<r separation of NH and IH might not be reached after one decade in the case of 
default systematic errors (Figs. 7 and 8), while it can be reached after ~ 6 years of data taking in the case of 
halved errors (Figs. 14 and 15). Similarly, assuming that the hierarchy is known, the energy-scale and flux-shape 
systematic uncertainties can significantly affect the accuracy of the (sf 2 , 5m 1 2 3 4 5 6 7 8 9 10 , ^m 2 e ) oscillation parameters emerging 
from prospective data fits (see Table I). 

The main message of our work is that further constraints on the admissible shapes and sizes of E —> E'(E) and 
4>(U) —> <&'{E) variations would be highly beneficial to the entire physics program of medium-baseline reactor projects. 
As side results of our analysis, we find the following: (1) the well-known IBD energy threshold acts as an effective self¬ 
calibration point in the fit; (2) neglecting matter effects may significantly bias the oscillation parameters (s 2 2 , 5m 2 ); 
(3) taking the wrong hierarchy may significantly bias the parameter Am 2 e ; and (4) prospective constraints on Th and 
U geoneutrino fluxes are largely insensitive to systematic uncertainties. 


Acknowledgments 

This work is supported by the Italian Istituto Nazionale di Fisica Nucleare (INFN) and Ministero dell’Istruzione, 
Universita e Ricerca (MIUR) through the “Theoretical Astroparticle Physics” projects. Preliminary results have 
been presented by A.M. at EPS-HEP 2015 , European Physical Society Conference on High Energy Physics (Vienna, 
Austria, July 2015). 


[1] K.A. Olive et al. (Particle Data Group), Chin. Phys. C 38, 090001 (2014). See the review therein: “Neutrino mass, mixing 
and oscillations,” by K. Nakamura and S.T. Petcov. 

[2] G. L. Fogli, E. Lisi and A. Palazzo, “Quasi energy independent solar neutrino transitions,” Phys. Rev. D 65, 073019 (2002) 
hep-ph/0105080 . 

[3] S. T. Petcov and M. Piai, “The LMA MSW solution of the solar neutrino problem, inverted neutrino mass hierarchy and 
reactor neutrino experiments,” Phys. Lett. B 533, 94 (2002) hep-ph/0112074 . 

[4] F. An et al, “Neutrino Physics with JUNO,” arXiv:1507.05613 [physics.ins-det]. 

[5] S. B. Kim, “New results from RENO and prospects with RENO-50,” arXiv:1412.2199 [hep-ex], to appear in the Proceedings 
of NOW 2014 , Neutrino Oscillation Workshop (Otranto, Italy, 2014), ed. by P. Bernardini, G.L. Fogli and E. Lisi, Nucl. 
Part. Phys. Proc. (Elsevier, 2015, in press). 

[6] A. Bandyopadhyay, S. Choubey and S. Goswami, “Exploring the sensitivity of current and future experiments to 
theta(solar),” Phys. Rev. D 67, 113011 (2003) hep-ph/0302243 ; S. Choubey, S. T. Petcov and M. Piai, “Precision neu¬ 
trino oscillation physics with an intermediate baseline reactor neutrino experiment,” Phys. Rev. D 68, 113006 (2003) 
hep-ph/0306017 ; A. Bandyopadhyay, S. Choubey, S. Goswami and S. T. Petcov, “High precision measurements of 
theta(solar) in solar and reactor neutrino experiments,” Phys. Rev. D 72 (2005) 033013 hep-ph/0410283 . 

[7] H. Minakata, H. Nunokawa, W. J. C. Teves and R. Zukanovich Funchal, “Reactor measurement of theta(12): Principles, 
accuracies and physics potentials,” Phys. Rev. D 71, 013005 (2005) hep-ph/0407326 . 

[8] P. Vogel, L. Wen and C. Zhang, “Neutrino Oscillation Studies with Reactors,” Nature Communications 6, 6935 (2015) 
arXiv: 1503.01059 [hep-ex]]. 

[9] X. Qian and P. Vogel, “Neutrino Mass Hierarchy,” Prog. Part. Nucl. Phys. 83, 1 (2015) arXiv: 1505.01891 [hep-ex]]. 

[10] S. H. Seo [RENO Collaboration], “New Results from RENO and The 5 MeV Excess,” Proceedings of Neutrino 2014, XXVI 
International Conference on Neutrino Physics and Astrophysics, ed by E. Kearns and G. Feldman, AIP Conf. Proc. 1666, 
080002 (2015) arXiv:1410.7987 [hep-ex]]! 



13 


[11] Y. Abe et al. [Double Chooz Collaboration], “Improved measurements of the neutrino mixing angle #13 with the Double 
Chooz detector,” JHEP 1410, 086 (2014) [Erratum ibidem 1502, 074 (2015)] arXiv: 1406.7763 [hep-ex]]. 

[12] L. Zhan [for the Daya Bay Collaboration], “Recent Results from Daya Bay,” arXiv: 1506.01149 [hep-ex], to appear in the 
Proceedings of NEUTEL 2015, XVI International Workshop on Neutrino Telescopes (Venice, Italy, 2015). 

[13] T. A. Mueller et al, “Improved Predictions of Reactor Antineutrino Spectra,” Pliys. Rev. C 83, 054615 (2011) 
arXiv:1101.2663 [hep-ex]]. 

[14] P. Huber, “On the determination of anti-neutrino spectra from nuclear reactors,” Phys. Rev. C 84, 024617 (2011) [Erratum 
ibidem 85, 029901 (2012)] arXiv: 1106.0687 [hep-ph]]. 

[15] D. A. Dwyer and T. J. Langford, “Spectral Structure of Electron Antineutrinos from Nuclear Reactors,” Phys. Rev. Lett. 
114, no. 1, 012502 (2015) arXiv:1407.1281 [nucl-ex]]. 

[16] A. C. Hayes, J. L. Friar, G. T. Garvey, D. Ibeling, G. Jungman, T. Kawano and R. W. Mills, “The Origin and Implications 
of the Shoulder in Reactor Neutrino Spectra,” arXiv: 1506.00583 [nucl-th]. 

[17] P. Huber, talk at the International Conference Neutrino Geoscience 2015 (Paris, France, 2015), available at the website 
www.ipgp.jussieu.fr/fr/evenements/neutrino-geoscience-2015-conference 

[18] F. P. An et al. [Daya Bay Collaboration], “A new measurement of antineutrino oscillation with the full detector configuration 
at Daya Bay,” arXiv: 1505.03456 [hep-ex]. 

[19] S. J. Parke, H. Minakata, H. Nunokawa and R. Z. Funchal, “Mass Hierarchy via Mossbauer and Reactor Neutrinos,” Nucl. 
Phys. Proc. Suppl. 188, 115 (2009) arXiv:0812.1879 [hep-ph]]. 

[20] X. Qian, D. A. Dwyer, R. D. McKeown, P. Vogel, W. Wang and C. Zhang, “Mass Hierarchy Resolution in Reactor Anti¬ 
neutrino Experiments: Parameter Degeneracies and Detector Energy Response,” Phys. Rev. D 87, no. 3, 033005 (2013) 
arXiv: 1208.1551 [physics, ins-det]]. 

[21] Y. F. Li, J. Cao, Y. Wang and L. Zhan, “Unambiguous Determination of the Neutrino Mass Hierarchy Using Reactor 
Neutrinos,” Phys. Rev. D 88 , 013008 (2013) arXiv: 1303.6733 [hep-ex]]. 

[22] F. Capozzi, E. Lisi and A. Marrone, “Neutrino mass hierarchy and electron neutrino oscillation parameters with one 
hundred thousand reactor events,” Phys. Rev. D 89, no. 1, 013001 (2014) arXiv: 1309.1638 [hep-ph]]. 

[23] F. Capozzi, E. Lisi and A. Marrone, “PINGU and the neutrino mass hierarchy: Statistical and systematic aspects,” Phys. 
Rev. D 91, 073011 (2015) arXiv: 1503.01999 [hep-ph]]. 

[24] A. de Gouvea, J. Jenkins and B. Kayser, “Neutrino mass hierarchy, vacuum oscillations, and vanishing - U(e3)—,” Phys. 
Rev. D 71, 113009 (2005) hep-ph/0503079 . 

[25] H. Nunokawa, S. J. Parke and R. Zukanovich Funchal, “Another possible way to determine the neutrino mass hierarchy,” 
Phys. Rev. D 72, 013009 (2005) hep-ph/0503283 . 

[26] H. Minakata, H. Nunokawa, S. J. Parke and R. Zukanovich Funchal, “Determination of the neutrino mass hierarchy via 
the phase of the disappearance oscillation probability with a monochromatic anti-electron-neutrino source,” Phys. Rev. D 
76, 053004 (2007) [Erratum-ibid. D 76, 079901 (2007)] hep-ph/0701151 . 

[27] F. Capozzi, G. L. Fogli, E. Lisi, A. Marrone, D. Montanino and A. Palazzo, “Status of three-neutrino oscillation parameters, 
circa 2013,” Phys. Rev. D 89, 093018 (2014) arXiv: 1312.2878 [hep-ph]]. 

[28] Y. Wang, talk at the International Conference on Massive Neutrinos (Nanyang Technological Univ., Singapore, 2015), 
available at the website: www.ntu.edu.sg/ias/upcomingevents/MassiveNeutrinos 

[29] D. V. Forero, M. Tortola and J. W. F. Valle, “Neutrino oscillations refitted,” Phys. Rev. D 90, no. 9, 093006 (2014) 
arXiv: 1405.7540 [hep-ph]]. 

[30] M. C. Gonzalez-Garcia, M. Maltoni and T. Schwetz, “Updated fit to three neutrino mixing: status of leptonic CP violation,” 
JHEP 1411, 052 (2014) arXiv: 1409.5439 [hep-ph]]. 

[31] K. Abe et al. [T2K Collaboration], “Neutrino oscillation physics potential of the T2K experiment,” Prog. Theor. Exp. 
Phys. 2015, no. 4, 043C01 (2015) arXiv: 1409.7469 [hep-ex]]. 

[32] V. Strati, M. Baldoncini, I. Callegari, F. Mantovani, W.F. McDonough, B. Ricci and G. Xhixha, Progress in Earth and 
Planetary Science 2, article ID 5, (2015) arXiv: 1412.3324 [physics.geo-ph]]. 

[33] D.A. Dwyer, talk at the Workshop “The Status of Reactor Antineutrino Flux Modeling” (Nantes, France, 2014), available 
at indico.cern.ch/event/353976 

[34] B.-Z. Hu, talk at Moriond EW 2015, 50th Rencontres de Moriond on ElectroWeak Interactions and Unified Theories (La 
Thuile, Italy, 2015), available at indico.in2p3.fr/event/10819 

[35] S. F. Ge, K. Hagiwara, N. Okamura and Y. Takaesu, “Determination of mass hierarchy with medium baseline reactor 
neutrino experiments,” JHEP 1305, 131 (2013) arXiv:1210.8141 [hep-ph]]. 

[36] In our opinion, it is appropriate to attach flux-shape uncertainties to the “theoretical spectrum” S and energy-scale 
uncertainties to the “experimental spectrum” S*. However, we have verified that our results are basically unchanged, if 
both uncertainties are assumed to act only on S or on S* . In such cases, in principle, one must also specify the ordering of 
the non-commutative operations E —>• E' and $ —> We have also verified that commuting such operations (on either S 
or S*) induces negligible numerical changes in our results. 

[37] M. Blennow, P. Coloma, P. Huber and T. Schwetz, “Quantifying the sensitivity of oscillation experiments to the neutrino 
mass ordering,” JHEP 1403, 028 (2014) arXiv:1311.1822 [hep-ph]]. 



